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A MATRIX METHOD FOR STUDYING MOTIONS OF SPACECRAFT 
CONSISTING OF INTERCONNECTED RIGID BODIES 

By Richard E. Turner 
Langley Research Center 

SUMMARY 

A system of matrix equations has been developed which describes the interaction 
between two rigid bodies joined at a point. The relative motion of the bodies is treated 
by projection matrices. A projection matrix is found which determines the angular accel- 
eration allowed by relative angular motion about a common attachment point. Another 
projection matrix is found which determines the components of relative angular acceler- 
ation of the two bodies caused by physical constraints on their relative angular motion. 

By using these matrices, an equation is found that relates the motions of the two bodies 
where physical rotational constraints are automatically applied. This equation, the 
translational momentum equation, the angular momentum equation, and the translational 
constraint equation form a set of equations which can be used to add rigid rotating bodies 
to a dynamical system to achieve a new dynamical system of any desired complexity. 

This method is useful when constraint forces and torques must be calculated. 

An example dynamical system which consists of a central hub with four rigid wings 
connected to the hub at four different attachment points is studied. The equations of 
motion for this system were integrated numerically, and some parameters which 
describe the motion of the system are presented graphically. 

INTRODUCTION 

Many spacecraft must erect or deploy in space panels such as solar cell panels and 
heat radiators or they must deploy booms for antennas or experiment mountings. The 
Pegasus type satellites, which deployed meteoroid penetration detectors with areas in 
excess of 200 meters^ are typical examples of these spacecraft. It is impossible to test 
these spacecraft in the zero gravity environment in which they must finally operate. Also 
since a spacecraft's physical geometry is usually determined by evolution from the orig- 
inal design to the final working design, it is too expensive and time consuming to build 
dynamic scale models for each evolutionary step. Thus, a need exists for an analytical 
method to calculate reliably the motions expected for such spacecraft while operating in 
space. 



Analytical studies of these large relative motion problems found in the literature 
usually fall into two groups. The first group contains those studies which use energy 
methods with emphasis being placed on the symmetry of the initial angular velocity of the 
dynamical system as well as on the physical symmetry of the system. References 1 
and 2 are examples of such studies. The results of these analyses are applicable only to 
dynamical systems having equivalent spin modes and physical geometry. The second 
group of studies are those that use vector mechanics. (See refs. 3, 4, and 5.) This group 
has results that are valid for general initial angular velocities. However, references 4 
and 5 are valid only for satellites which deploy simple booms (booms consisting of single 
lumped mass points). The results of reference 3 are valid only when the motion of the 
central body is known. Thus neither group of studies gives results which are generally 
valid for the calculation of motions for satellites of a general geometry. 

A method for calculating the motions of spacecraft consisting of interconnecting 
rigid bodies is presented in this report. The method can be applied to structures of any 
general configuration or changing configuration such as a spacecraft that deploys com- 
ponents. The method allows the motions of the rigid bodies to be calculated as well as 
the interaction forces and torques between the bodies. The number of bodies in a single 
system which can be studied by this method is limited only by the speed and storage of 
modern digital computers. Matrix notation is used because it is thereby possible to cast 
the governing equations in terms of variables common to all such dynamic systems. The 
resulting equations can therefore be used (as is) to study any such dynamical system. 

An application of the method is presented in which the motion of a proposed NASA 
micrometeoroid spacecraft deploying four large panels is studied. Some motion time 
histories are generated for this dynamic system which consists of five rigid bodies - 
the central spacecraft body plus four rigid bodies representing the deploying panels. 

Two means of deploying the panels are investigated. In the first case the panels are 
assumed to be deployed by torsion bars with viscous damping. Motions are calculated 
for situations for two damping rates and for situations in which one panel is allowed to 
deploy before the other three. The second case studied assumed that the panels are 
symmetrically deployed by centrifugal force resulting from spin of the spacecraft. This 
second case serves as a check on the method since the angular momentum of the five- 
body system is conserved. The analysis given here is valid up to panel locking only, but 
it is possible to extend the method to cover such cases as panel locking. 

SYMBOLS 

The notation shown is similar to the Dirac notation used in reference 6. 
a acceleration of center of mass of body, meters/second^ 
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A 

A(i,j) 

B 

B m 

C 

C Y 

D(i) 

e 

E(i,j) 

f 

F 

G 

H 

U 

H 

L 

L a 

L c 

m 


moment of inertia about center of mass of body, kilogram-meters^ 

a set of equations governing motion of two interconnected rigid bodies 
(body i and body j) when no rotational constraints are present 

moment of inertia about reference point of a body which does not coincide 
with its center of mass, kilogram-meters^ 

similar to B and nonsingular (see eq. (27)), kilogram-meters^ 

coordinate frame 

inertial coordinate frame 

a set of equations governing motion of body i 

a unit base for an orthogonal coordinate frame, dimensionless 

a set of equations governing motion of body j relative to body i 

net force on a body, newtons 

interaction constraint force, newtons 

force on a body, newtons 

component of angular momentum, kilogram-meter s^/second 
denotes specific bodies 
identity matrix 

interaction torque, newton-meters 
applied interaction torque, newton-meters 
constraint interaction torque, newton-meters 
mass of a body, kilograms 
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M 

N 

Pa 

Pc 

q 

Rl,R2,R3 


S(1,2,...,N) 


t 

T 

V^V 2 

V(P) 

W 

X 

Y,Z 

[»> 


a 


(3 





transformation matrix, dimensionless 

torque on a body, newton-meters; also integer 

projection matrix for applied interaction torque, dimensionless - 

projection matrix for constraint interaction torque, dimensionless 

inertial acceleration of reference point on a body, meter s/second^ 

rotation matrices (see eqs. (29) to (31)) 

a set of equations governing the motion of an N body dynamical system 
which is composed of interconnected rigid bodies 

time, sec 

total torque applied to a body about its center of mass, newton-meters 
vector in and C 2 representation 

velocity of reference point P in inertial space, meters/second 
damping rate 

matrix made up of components of three vectors 
elastic torque rate constants, newton-meters/radian 
column vector with all components being zero 

vector from a body’s reference point to its center of mass, meters 
vector between reference points of two bodies, meters 
Kronecker delta 

modified Euler rotation angles (see fig. 5), radians 
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V 


unit vector 


c«j angular velocity, radians/second 

S2 position coinciding with a, meters 

£ a vector which is reciprocal to rj, dimensionless 


Notation: 


<] 

row vector 

[> 

column vector 

PI 

length of £ y 

[] 

square matrix 

LI 

antisymmetric square matrix 

□ T 

transpose of [ J 

[]-* 

inverse of £ J 

(-) X (-) 

vector cross product 


Subscripts designate the body to which a given property belongs and superscripts 
designate the coordinate frame in which the property is represented. Bars over symbols 
denote vectors. A dot over a symbol denotes a time derivative with respect to the coor- 
dinate frame denoted by the superscript and ^ denotes a time derivative with respect 

dt 

to an inertial coordinate frame. 

ANALYSIS 

This method is not applicable to ideal locking joints (joints where the torque rates 
instantaneously change from finite to infinite values) during locking time. However, ideal 
locking joints may be approximated satisfactorily with highly damped, high-spring-rate 
joints during locking periods. After joint locking occurs, this method is applicable with- 
out approximation. 
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Mathematical Model 


The dynamical system model under consideration is one which has gross changes 
in its physical geometry and one which can be constructed by adding a single rigid body 
with a single attachment point to form a new system from an old one. An example of 
such a dynamical system is a spacecraft with large folding panels. The motions of such 
a system can be analyzed by considering the motions of the bodies about the attachment 
points. When a dynamical system consists of only one body, that system's motion is 
described by the path of the center of its mass and its rotational motion. If a second 
body is added to this system at some common attachment point, the interaction of the 
second body with the first body must be found. There will be a force vector and a torque 
vector present which reflect physically the attachment of the two bodies. When a third 
body is added to the system, its interaction with one of the first two bodies must be found. 
Again there will be force and torque vectors that show how the third body is added to the 
older system to form a newer system. By this procedure a dynamical system of N 
rigid bodies can be constructed and studied. 

Outline and Discussion of Method 

Dynamical systems are governed by second-order differential equations where 
position coordinates and velocity coordinates must be specified before the accelerations 
can be calculated. Although the governing equations of motion may be nonlinear functions 
of position and velocity coordinates, these equations are always linear functions of the 
accelerations. In general, for free dynamical systems of N rigid bodies there are 
N vector angular acceleration and N vector translational acceleration equations. For 
the dynamical systems considered here, there are always N - 1 vector translational 
constraint equations which can be used to reduce the number of vector equations from 2N 
to N + 1 with N + 1 vector unknowns, these being N vector angular accelerations and 
one vector translational acceleration. Rotational constraints may also be present in the 
system, and thus reduce still further the number of equations to be solved simultaneously. 
When the number of equations has been reduced to the smallest number, the mathematical 
similarity between different systems is lost. The method of solution given here uses a 
different approach. The advantages of this method will depend upon the dynamical system 
involved. Here the emphasis is placed on the similarity between systems and the analy- 
sis is carried as far as possible without using detailed knowledge of the system contraints. 

The method of solution defines the motion of an arbitrary body in the system in 
terms of that body's interaction with the system. The body's motion is described by its 
translational and angular acceleration, whereas its interaction with the system is 
described by the forces and torques impressed on the system by the body. Because all 
these vector equations are linear in accelerations, forces, and torques, these quantities 
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can be found by classical matrix algebra. The matrix notation used here is introduced 
in appendix A. 

The actual solution for a given system can be obtained by forming an independent 
set of vector equations equal in number to the total number of translational and angular 
accelerations as well as the interaction forces and torques, and then solving this set of 
vector equations for the unknowns. The angular accelerations are then integrated to find 
the angular coordinates and velocities at a new time. The translational accelerations 
need not be integrated because the coefficients in the governing set of equations are func- 
tions of angular variables. 


Solution of One and Two Body Systems 
With No Rotational Constraints 


The equations of motion which govern the motion of one and two body systems are 
developed in appendix B. Let Cj represent the coordinate frame fixed in body i, where 
body i is a single body system. The unit orthogonal base vectors for are e^(l), 
ij(2), and e^(3). From appendix B the translational momentum equation for body i is 


mi 



( 1 ) 


where the subscripts designate the body to which a given property belongs and super- 
scripts designate the coordinate frame in which the property is represented and mj is 
the mass of body i, £a) is acceleration of the center of mass and [f ) is the total 
external force on the body. The brackets £ ) denotes column vectors. The angular 
momentum equation for body i taken from appendix B is 


[4]^Y>=[Ti)-L"M[#iY> » 

where is the angular velocity and the double subscript indicates that 

rotation of relative to Cy where Cy is the inertial coordinate frame. [A] is 
the moment of inertia matrix about the body’s center of mass, |_wj is the antisymmetric 
matrix made up from [w) analogous to the w x vector operator as shown in appen- 
dix A and [t) is the total external torque on the body about its center of mass. These 
are the equations of motion for a single rigid body system. 

Consider a system composed of two rigid bodies which are attached at a common 
attachment point (reference point of body j) as shown in figure 1. Equations (1) and (2) 
are the equations of motion for body i. The translational momentum equation for body j 
taken from appendix B is 


L-5575 
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( 3 ) 



where jjq^ is the inertial acceleration of the reference point of body j, is the 

vector from the reference point of body j to the center of mass of body j. The angular 
momentum equation for body j, measured about its reference point, taken from appendix B 
is 







(4) 


where [bQ is the moment -of -inertia matrix measured about the reference point of the 
body and (l) is the total torque on the body about its reference point. The translational 
constraint equation for the attachment of body j to body i taken from appendix B / and 


noticing that 



equals 




(5) 


where [m] is a transformation matrix and its double subscript indicates that it con- 
verts column vectors in Cj representation to Cj representation. is a vector 

from the reference point of body i to the reference point of body j. In every case a ref- 
erence point of a body is the same as its attachment point if it has an attachment point. 


Next consider the interaction force and torque at the reference point of body j and 
designate as [F(ij)^ and [L(ij)^, respectively. The term Jljj^ from equation (4) can 
be rewritten as 


and 


Lj)= [Lj(i^ ♦ [ N j) 

f from equation (3) can be rewritten as 



( 6 ) 


(7) 


In equation (6), is the external torque applied to body j about its reference as 

opposed to the interaction torque jLj(ij)^, and in equation (7), jG^ is the external 

force applied to body i and |Fj(ij)^ is the interaction force. The term [ T i) in equa 
tion (2) can likewise be rewritten as 
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( 8 ) 



and 



in equation (1) can be rewritten as 



( 9 ) 


Equations (6) to (9) can be substituted into equations (1) to (5) to obtain the following 
set of equations A(i,j): 



where equation (14) was obtained from equation (5) by multiplying equation (5) by [^jjj • 
Equations (10) to (14) are the equations of motion for a dynamical system composed of 
two attached rigid bodies where no rotational constraints exist. The terms on the right- 
hand side of these equations must be known. 


For the case of a single body system 


, Fj(ijj) 


and 


Lj(ii) 


are zero and thus 


equations (10) and (11) remain as the describing equations where (jST^ and [g^ are 
prescribed forcing functions. Then and J^yy> can be calculated. 

For the case of a two-body system, equations (10) to (14) are the five governing 
vector equations. The vector unknowns are Jaj^), fe)- [“Sy) “ d Wefl). 

It can be shown that this set of equations does have a unique solution for the unknowns. 
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Application of Rotational Constraints 

Consider the case in which the angular velocity of body j relative to body i is con- 


strained to be zero about one or two axes and 


44 


can no longer be specified as a 


function of position and angular velocity. The equation set A(i,j) defined by equations (10) 


to (14) can no longer give the required solution because has now become an 


unknown. Let 


Laj ( 4 



be the known interaction torque specified by position and angu- 


lar velocity; also let 
Then 


h<4 


be the unknown constraint torque of body j on body i. 


|bij (ij^> = Laj(ijj) + |jjCj(ij)) 


(15) 


Two cases are considered. In the first case, body j is allowed to rotate relative to body i 
about two axes; in the second case, body j is allowed to rotate relative to body i about 
one axis. 


Relative rotation about two axes .- In the first case consider three unit vectors 
^1(1^ [^j(2^> 311(1 jr/j(3^ such that the first two form the instantaneous axes (not 

necessarily orthogonal) of relative rotation of body j to body i and the last one is 



The component Lj(ij^ is known along ^(1^ 311(1 7?j(2^. The unknown component 


Le](iij> 


must therefore lie along 


~ 4 3 >) 


and must be perpendicular to 


4 ») 


and 


jr?j(2^>. Consider also the vectors |?:j(l^>, j?j(2^, and j?j(3^> which are reciprocal 

to [r?j(lJ>, |^j(2^, and [^(3)) as defined in appendix C. Let jp a jj 311(1 j^c jj be 


defined by 



( 17 ) 
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and 


H ’ B (3 K ? i (3) ] <18) 

The matrix |^? a jj is the projection matrix for the angular velocity of body j relative to 

body i, and for |l. a |(ij^. Projection matrices and reciprocal vectors are 

defined in appendix C. The matrix jp c jj is the projection matrix for |lj C j(ij^. The 
vector j^°ji^ ln c j representation is given by 

Ri) = Ry) - [ M ij] Ry) 

It follows from equations (17) and (18) that 

[Paj] [l]«)> = [Laj(U)> 

and 

N][4>=[°> 


(19) 


( 20 ) 


( 21 ) 


Equations (20) and (21) are the additional equations needed to solve the problem where 
the rotational motion of body j relative to body i is constrained to lie along two axes. 


Relative rotation about one axis.- Let the axis of relative rotation be 
this case 


3, Lj(ijj) along r?j(l)^> i 



;i 

and 


Because 


L c ](ij)) i 


is known so that 
is a linear combination of 


are both perpendicular to 


to jr?j(2^ and £?j(3)), j??j(lj) 311(1 

M is 


r?j(lj). Si 


Since 


[r?]( lj). In 

is perpendicular to 
7 ]j( 2 ^> and |^?j(3^>, |j?j( 2 j) 
?|(1^ is also perpendicular 


must be parallel and equal. Thus, 


is defined by 



( 22 ) 
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and since the sum of 
is always defined by 



and 



is equal to the identity matrix 

[I] - M 



and |l? a jj and |^P c jJ are again known * Equations (20) and (21) can now be used to 
impose the desired constraints on rotational motion of body j relative to body i. 


Return to equation set A(i,j) and recall that when j^Lcj(ij)) is different from [0^, 

equation set A(i,j) has six vector unknowns and only five vector equations. One more vec- 
tor equation or three scalar equations are required to obtain values for the unknowns. If 
body j is allowed to rotate relative to body i about one (or two) axes, equation (20) pro- 
vides one (or two) new independent scalar equations. Also the time derivative of equa- 


tion (21) provides two (or one) new independent scalar relations between and 

jd)|^. Therefore equation (20) plus the time derivative of equation (21) provide the three 

necessary new scalar equations necessary to solve equation set A(i,j) for the unknowns. 
These three independent scalar relations may be put into vector form simply by adding 

the time derivative of equation (21) ^ premultiplied by to obtain the correct order 

of magnitude^ directly to equation (20) to get 



Equation set A(i,j) plus equation (24) can now be solved directly by matrix inversion to 
get unique values for the unknowns. However, if one desires to find the unknowns by 
algebraic substitution (this will generally be the case because direct matrix inversion 
will require more computer storage and longer computation time than solution by alge- 
braic substitution), a useful alternate to equation (24) can be found by premultiplying 

equation (13) by jp a jl and adding the result to equation (24) then 
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( 25 ) 



which can be shown to be nonsingular. When equation (26) is substituted into equation (25), 
there follows 



Since L B m ij is nonsingular, equation (27) 
bers of set A(i,j) to get unique solutions for 



can be combined algebraically with the mem- 



Extension to multibody systems .- Consider now an N body dynamical system. 
From these bodies, arbitrarily select a reference body, body i, and let its reference 
point be its center of mass. The motion of this body is governed by equation set D(i) 
defined as equations (10) and (11). Any other body, body j, which is attached to the ref- 
erence body has its attachment point as its reference point and its motion is governed by 
equation set E(i,j) defined as equations (12) to (14) plus equation (27). Any other body, 
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body k, which is attached to a body j has its attachment point as its reference point and 
its motion is governed by equation set E(j,k) and so on. Thus the reference body, body i, 
has its motion governed by D(i) and every other body, body m, has its motion governed 
by E(Z,m) where body l is the body to which body m is attached and the reference point 
of body m is the point at which it is attached to body l. The set of equations which gov- 
erns the motion of the entire system S(1,2,...,N) is D(i) plus all the E(i,m) for 
every other body m in the system. When the physical properties of the system along 
with initial coordinates and coordinate velocities have been specified, the system’s 
unknowns can be found by classical matrix inversion or by algebraic substitution between 
the members of S(1,2,...,N) according to the desires of the investigator. Any desired 
number of unknowns may be calculated without calculating the remaining unknowns. If 
one is to generate a time history of the motion of this type of dynamical system, at least 
all the angular accelerations must be calculated to provide angular positions and angular 
velocities which in turn determine the coefficients in the equation set S(1,2,...,N). A 
five-body system is illustrated in the example problem. 

One word of caution is in order. When the system's motions are solved from 
S(1,2,...,N), the correct accelerations will automatically be found. However, when these 
accelerations are integrated, one must be exceedingly careful to insure that the integra- 
tion scheme conforms to the motion constraints applied to the system. This precaution 
is necessary because all numerical integrations schemes have some inherent error and 
over long time intervals, motion constraints might otherwise become grossly violated. 

APPLICATION 

Translation of a Spacecraft's Physical Properties into Matrix Variables 

The analysis is illustrated by studying the deployment motions of a proposed NASA 
micrometeorite detection satellite shown in figure 2. This satellite would essentially be 
composed of five rigid bodies connected at four attachment points. The body which would 
house the satellite's communication electronics (the hub) is centrally located between 
four other bodies (the wings) which are the meteorite detection sensor arrays. Each 
wing is made up from three smaller rigid members, a main panel plus two auxiliary 
panels. In the launch configuration the auxiliary panels are folded down to the main 
panel and all four main panels are folded about the ends adjacent to the hub into a square 
cylinder, the hub forming a cap for the cylinder, as shown in figure 2(a). 
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The initial phase of deployment consists of the four wings folding out of the cylin- 
drical configuration (with the auxiliary panels still locked to the main panels) into a 
planar configuration as shown in figure 2(b). Then the auxiliary panels are unfolded into 
the position shown in figure 2(c). 

The motion studied here encompasses the time from which the wings begin to 
deploy from the cylindrical configuration until they are in the planar configuration. 

Since the auxiliary panels (during this period) are folded against the main panels and 
locked in that position, the wings are assumed to be rigid bodies, and the whole system 
is to be treated as the ideal system shown in figure 3. Each wing is assumed to rotate 
relative to the hub about two axes. (See figs. 3 and 4.) The first axis of rotation, the 
hinge line, is the axis of wing deployment. The second axis of rotation is perpendicular 
to the hinge line and passes through the center of mass of the individual wing. This 
degree of freedom is intended to provide a pseudo-simulation of any twist which develops 
in the wings during deployment. The coordinate frame for the hub is shown in figure 4 
along with the hinge line for body 2 and the line through the center of mass of body 2 
about which body 2 can rotate relative to body 1. 

Figure 5 gives the modified Euler angles which relate (coordinate frame fixed 

in body 1) to C 2 (coordinate frame fixed in body 2). Consider an intermediate coordi- 
nate frame C which initially is alined with Cj. First, rotate about e(3) counterclock- 
wise through angle i/^. Second, rotate about the new e(l) counterclockwise through 
angle 62 . Finally, rotate about the new e(2) counterclockwise through angle 72 • Tiie 
intermediate coordinate frame is now alined with C 2 . It follows that a vector in 

representation [v^ is related to the same vector in C 2 representation [v^) by 

fr 1 ) = [ R l] [ R 2 ] [ R 3 ] [I 2 ) - [ M 2l] I? 2 ) < 28 > 

where 


[ r j = 


cos 1^2 
sin 1^2 
0 


-sin 1^2 
cos 

0 


0 

0 

1 


(29) 
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1 0 

0 

p*]“ 

0 cos 02 

-sin 0 2 


0 sin &2 

cos 0 2 


cos 7 2 0 

sin 7 2 

II 

l£Ll 

0 1 

0 


-sin 7 2 0 

COS 7 2 


( 30 ) 


(31) 


lotation may occur about two axes, ^(l) and rj 2 (2)- Figure 5 shows that 77 2 (1) (the 
xis of 0 2 rotation), 7? 2 (2) (the axis of y 2 rotation) and fyg (3) (the axis along which 
lie constraint torque lies) form a set of mutually orthogonal unit vectors. In all such 
ases, the vectors 5 2 (1) (which are reciprocal to the vectors 772 (i)) are equal to the 

ectors 772 ( 1 )- From figure (5) it can be seen that 


772(1) = e 2 (l) cos y 2 + e 2 ( 3 ) sin y 2 
r? 2 ( 2 ) = e 2 ( 2 ) \ 

77 2 ( 3 ) = -©2(1) sin r 2 + e 2 ( 3 ) cos y 2 


(32) 


t now follows directly from equation (17) that 


1.6 




(cos r 2 ) 2 0 

0 1 

COS sin y 2 0 



cos 7 2 sin 7 g 
0 

(sin 7 2 ) 2 


< 0 , 1 , 0 ] 


( 33 ) 



also since 



Then from equation (18) 



(sin y 2 ) 2 
0 

-cos y 2 sin y 2 


0 -sin y 2 cos y 2 
0 0 

2 

0 (cos y 2 ) 


(34) 


The three other wings (bodies 3, 4, and 5) are physically connected to the hub in the same 

for each wing can be 

found for each wing by replacing ^yg, ^ ^2) by ^j) * n e Q. ua ^ ions (29), (30), (31), 

(33), and (34). The angles i/a give the spacing of the jth wing around the hub, whereas 
0 j and y^ represent the degrees of freedom of each wing relative to the hub. The 

torques applied to each wing were assumed to be applied at the point of attachment for the 

attachment of the wing to the hub. One torque component was assumed proportional to both 

0. and 0. and parallel to the axis of rotation for 0.. A second torque component was 
3 3 3 

assumed to be proportional to y. and parallel to the axis of rotation of y-. The physi- 

j j 


manner as body 2; thus, the matrices 


p J p ] 

_ a 3j ’ L C LT 


and 


cal properties of the example spacecraft are given in table I. The vector, in Cj 

is measured from the center of mass of the hub to the attachment point of the jth wing. 
The vector o>^) is measured in Cj coordinates from the attachment point of the hub 
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TABLE I.- PHYSICAL PROPERTIES OF EXAMPLE SPACECRAFT 


Radius vector from attachment to center of mass for wings 2, 3, 4, and 5, dimension, m: 








Radius vector from reference point of hub to attachment point, dimension, m: 



Moment-of -inertia matrix for hub about its center of mass, dimension, kg-m2: 

877 0 0 “ 

0 877 0 

0 0 1755 



Moment-of -inertia matrix for jth wing about its attachment point to hub, 
dimension, kg -m2 



83000 

0 

0 


0 0 

2610 0 

0 85700 
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to the wing to the center of mass of the wing. The mass of the hub and wings were each 


2620 kg. The moment of inertia matrix for the hub about its center of mass 



is 


given in table 
points |b| 


I. The moment of inertia matrices for all wings about their attachment 
were identical to each other. All four wings are identical and their respec- 


tive coordinate frames are each oriented relative to each body in identical arrangements 
and this arrangement is given in figure 6. The spacing of the wings around the hub is 
given by i^ 2 , ^3, 1^4, and equal to 0, tt/ 2, it, and 3ir/2 radians, respectively. The 

torques applied to the wings about their point of attachment are governed by the torque 
coefficients Za, W i? and Y i for each wing: Z, is the torque rate (N-m-rad"l) about 

the axis of 0j rotation, Wj is the torque damping rate (N-m-sec-rad~l) about the axis 


of (9j rotation, and is the torque rate (N-m-rad - l) about the axis of yj rotation. 
Thus, jL a |(lj^ (the torque applied to body j at its point of attachment) is given (in 
Cj coordinates) by 



(35) 


When the initial spin rates 


“jy) 


of all the bodies in the system are specified 


along with 1^,, 0., and y. for each of the four wings, equation (28) is used to calculate 

J J J — V 

[ M »] for all four wings, after which equation (19) is used to calculate 


for each 


wing. From figure 5 it is seen that 



(36) 


Since 




ovL\ is known for each wing, 0. and y , can be found for each wing and 

_ J J J 


can be calculated by using equation (35). Forces and torques are applied to 


the wings and hub only at the attachment points so that since the external forces and 
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torques are all zero, |g g|^), and [ N l) in equations (10) 

zero. Thus the equation set (D(l)) o-nvprnino- thp mritinn nf thp huh is 


to (13) are all 


r 


M 


1 



Set D(l) = J 


j=2 



v. 



+ 



( 37 ) 



and the equation set E(l,j) is 



(42) 


The equation set S(l,2,3,4,5) which governs this dynamical system is given by 

S(l,2,3,4,5) = D(l) + E(l,2) + E(l,3) + E(l,4) + E(l,5) (43) 


This set of simultaneous equations can be solved by classical matrix inversion or, since 
all the equations in the set are three dimensional, they may be solved by algebraic sub- 
stitution within the set. 

The system of equations defined by S(l,2,3,4,5) in equation (43) was programed 
for a digital computer. This linear system of equations was solved by algebraic 
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substitution to get unique values for j^a^, [qj),and 


Then tejw and dij^) 


were used to determine 


JwJ Y ^) and jw| Y ^ at 


new time, and equation (36) was used to update 8^ and y^ so that [“.»] could be 
updated. Equation (A4) was used to update j^iy] in the form of iteration by 


' E MlY Jt * 5 "([ m iy] + ([“«] 


At 


(44) 


and it was orthogonalized at each time increment by averaging [ m iy] with 

(W l ) • This cycle was repeated until the desired time history of the system was 
produced. A flow chart for the computation sequence is given in table n. 

It should be noticed that to appears as a quantity which connects the hub back 
to its initial coordinates in inertial space and it need not be retained if one desires only 
the time histories of all jMj^J . 

The results of three computer runs are presented here. These runs correspond to 
two different methods of deploying the spacecraft's wings from the launch configuration. 


Motions From Nonsymmetrical Deployment With No Spin 

One method of deploying the wings is initially to have the spacecraft nonrotating in 
space and use torsion bars with viscous damping to fold the wings into the planar config- 
uration. This method of deployment was studied for two damping rates to establish the 
effects of damping on nonsymmetrical deployment motions. 

Case 1 is the simulation of this type of deployment with small damping in the 

deployment system. This case shown in figure 7 has the initial conditions Hr) 
equal [O^ for all four wings and hub, and 0j equal ir/2. The torque rate about the 
axis of rotation Zj for each wing was 482 N-m-rad - l, the torque damping rate 
about the axis of 0j rotation Wj for each wing was 560 N-m-sec-rad~l, and yj was 
constrained to be zero. Wing 2 was released 2 seconds before wings 3, 4, and 5. The 
time integration was continued until wing 2 rotated from 82 equal n/2 to 82 equal 
zero. The integration was stopped at that time to avoid the event of panel locking. Fig- 
ure 7 (time histories of 82, 83, 84, and 05) shows that the 2-second delay caused an 

angular lag for the wings released late. 
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TABLE H.- FLOW CHART OF COMPUTATION SEQUENCE 


1. Input physical properties: 


[4 K 

[H> K> 

(j = 2,3, 4, 5; i = 1) 

[V l) > 

(j = 2, 3,4, 5; i 

= 1,2,3) 

m., : 

b w i- Y i 

(j = 2, 3,4, 5) 

2. Input initial conditions: 

[ 

m 1y ]; [»ir) 

(i= 1,2, 3, 4, 5) 

8j, Yy 

0 = 2, 3,4, 5) 


3. Calculate all coefficients and constants in the linear set of equations S(l,2,3,4,5) at 

the initial time. 

4. Calculate all the unknowns in the set of linear equations S(l,2,3,4,5) by algebraic 

substitution at the initial time. 


5. On first iteration cycle, set the acceleration at the new time equal to the accelerations 
at the initial time. 


6. Integrate accelerations to get the position coordinates and velocity coordinates at the 

new time by the formulas: 

% + it=«t + «*t at+ |(r + V s ) A * 2 

^t+At = Of + (§t + ^t+At)f^ 

7. Calculate all coefficients and constants in the linear set of equations S(l,2,3,4,5) at 

the new time. 

8. Calculate all the unknowns in the set of linear equations S(l,2,3,4,5) by algebraic 

substitution for the new time. 

9. Cycle N times between points 6 and 8. 

10. Output desired quantities; replace initial conditions with new time position coordinates 
and velocity coordinates; go to point 5. 
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Case 2 is the simulation of this type of deployment with large damping in the system. 
For this case the initial conditions and constants were the same as those of case 1 except 
that the damping torque rate about the axis of 0j rotation Wj for each wing was 
5530 N-m-sec-rad"l. Figure 8 (time histories of &2> &3> #4> an( * 65 ) shows again 

that there is an angular lag for the wings released late and the deployment time is 
increased about 50 percent. 

Figure 9 gives a comparison of the pitching motion of the hub induced by the non- 
symmetrical deployment of the wings for cases 1 and 2. The angle time histories plotted 
correspond to the angle between the vector e^(3) at time zero and the vector ej(3) 
for times greater than zero, where e^(3) is the third axis of the coordinate frame 
fixed in body 1. Figure 9 shows that increased damping (in the deployment system) 
decreased the pitching motion of the hub when a nonsymmetrical deployment occurs. 

This result probably would not hold true for dynamical systems different from the one 
studied here. 

In order to establish a check on the digital computer program used to calculate 
time histories of the example problems, an independent exact solution was calculated for 
a symmetrical deployment in which no damping was present and the results are plotted 
in figure 10. Figure 10 is in good agreement with figure 7. It should be remembered 
that both runs had no spin but the solution given in figure 7 had a small amount of 
damping. 


Motions From Symmetrical Deployment With Spin 

Another method of deploying the wings is to use a combination of damped torsion 
bars with vehicle spin about ej(3) so that centrifugal force on the wings speeds up the 
deployment process. Case 3 shown in figures 11, 12, and 13 gives such a deployment 

time history. This case has the initial conditions e< l ua l * or ^ ^ our w i n S s 

and Ry> equals 10 revolutions per minute about e^(3). Initially, the deployment 
angle 0j was n/2 and the twist angle was zero for all four wings. For this case, 
y- was not constrained to be zero; the elastic torque rates were each 

4330 N-m-rad“l; the elastic torque rates Zj were each 482 n-m-rad~l; and the damping 
rates Wj were all zero. Figure 11 gives the time history of the deployment angle 0j 

from which it can be seen that the deployment time has been appreciably shortened rela- 
tive to the no spin case of figure 7. Figure 12 gives the time history of the wing twisting 
which develops as the wings deploy. The size of the twist angle suggests that a stiffer 

wing is needed or else a smaller initial value of Ry) should be used. Figure 13 
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^gives the third component of 



and demonstrates the manner in which the space- 


craft spin slows down as the wings deploy. 


This case also serves as a numerical check on the correctness of the overall 
approach because it is easy to verify that the angular momentum of the system about 
ei(3) has been conserved during the deployment of the wings as it should. In all cases, 
i//j for all four wings was constrained to remain constant. 


CONCLUDING REMARKS 


The motions of dynamic systems built up of interconnected rigid bodies, which 
undergo large relative angular displacements, can be calculated with matrix calculus. 
The method given in this report makes no use of such quantities as "generalized 
coordinates" and therefore applies (as is) to any dynamical system of interconnected 
rigid bodies. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., February 20, 1968, 

124-09-14-04-23. 
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APPENDIX A 


PROPERTIES OF SPECIAL MATRICES 

Consider an orthogonal coordinate frame, designated C, and use subscripts to 

refer to specific coordinate frames. If a vector's components [v‘> are given in Cj, 
the same vector's components in Cj are given by reference 7 as 

[v>)= [My] [v 1 ) (Al) 

For the double subscript shows the direction of the change in the coordinate 

frames; in this case the subscript ij shows that the C^ representation is being 
changed to the Cj representation. 

Next consider two vectors A and B and their cross product C so that 

A x B = C 


If 



and if an antisymmetric square matrix j_A~( is formed from the components of [a), 
so that 






it is simple to verify that 



APPENDIX A 


The scalar product of 


[A) and 


is given by 


<a][b> = (b][a> 


(A3) 


where (aJ is the transpose of [a). Finally consider again coordinate systems 
and Cj and the orthogonal transformation matrix j^ijj • The time derivative of 

K] is given by (see ref. 7) 



(A4) 


where 


as in equation (A2) and 


[4] is the antisymmetric matrix made from j^ij^ 

[4) is the angular velocity vector of relative to Cj and represented in C^. 
For orthogonal matrices, the inverse is equal to the transpose; thus, 


N' 1 = [ M iif <As 

The matrix notation used here is similar to Dirac notation (see ref. 6) and it has 
the property that the outermost brackets in the product of two matrices always indicates 
the matrix character of the product. 
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APPENDIX B 


THE TRANSLATION OF VECTOR EQUATIONS OF MOTION 
OF RIGID BODIES INTO MATRIX NOTATION 

The equations of motion for a rotating rigid body are 


ma = f 


(Bl) 


and 


d(H) = f (B2) 

dt 

where m is the body mass, a is acceleration of the center of mass, f is total force 
applied to the body, H is total angular momentum about a reference point, L is the 
total torque applied at the same reference point, bars over quantities denote vector char- 
acter and denotes time differentiation of the parenthetical quantity in the inertial 

dt 

coordinate frame. 

Consider body i in figure 1 and let its reference point be its center of mass. Equa- 
tion (1) then becomes 


m i^i = * i ( B3 ) 

where the subscripts denote that a property belongs to body i and superscripts denote 
the coordinate frame in which the nonscalar quantities are measured. Equation (B2) 
becomes 


n\ + wj Y x Hj = Lj (B4) 

where H denotes time differentiation of H in a body -fixed coordinate frame (denoted 
by the superscript) and lo is angular velocity and its double subscript indicates that it 
measures the spin of coordinate frame i relative to the inertial coordinate frame Cy 

Next consider body j in figure 1 and let its reference point be point P (the point of 

attachment of body j to body i). In this case aj (the acceleration of the center of mass 

of body j) can be written as 

5 M-s-( 5 ! x “!y) (b5) 
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APPENDIX B 


where q. is the acceleration of the reference point of body j and a?, is the position 
] ** 
vector measured from the reference point of body j (fixed in body j) to the center of mass 

.of body j. When equation (B5) is expanded, 


aj = qj - wj Y x (ai x - aj x wj (B6) 

When equation (Bl) is written for body j in terms of its reference-point acceleration, the 
result is 


m i [ 5 i ' a i x - 4 x ( 5 S x “jy) 


- f S 


(B7) 


The angular momentum equation for "body j" about its reference point is 


dt\ J 


±m) = £i 


(B8) 


Equations (B3), (B4), and (B7) can be put in matrix notation by replacing all vectors and 


vector operations by matrices as in appendix A. In equation (B4) let h| be replaced 


by [Hj) 


and then let 


GO 


be replaced by 


a[| [ 4 ) 


where 


is the inertia 


matrix of body i about its center of mass. The results are 



and for equation (B7) 


m. 




Equation (B8) can be transformed into matrix notation after noticing that 


(B9) 


(BIO) 


(Bll) 




“j 




(B12) 


where JjBjj is the inertia matrix of body j about its reference point and y is a 
position vector measured from a point in inertial space (coinciding with its reference 
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point) to the center of mass of body j. It must be noted that although S2. and a. are 

3 J 

equal at the time instant under consideration, their derivatives are not equal and thus a 
distinction must be made whenever differentiation is performed. When equation (B12) 
is substituted into the matrix equivalent of equation (B8) and after performing the indi- 



Finally consider the constraint caused by body j being attached to body i. The 
resulting equation in terms of accelerations is 

* i + f ( s 1y x < B14 > 

where js. is the vector from the reference point of body i to the reference point of 
body j as shown in figure 1. When the indicated differentiation in equation (B14) is per- 
formed, there follows 



which, when translated into matrix notation, gives 


(B15) 



where * s the square matrix that changes column vectors represented in the 

jth coordinate system to the same column vector represented in the ith coordinate system 
as shown in appendix A. 
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APPENDIX C 


PROPERTIES OF PROJECTION MATRICES 


Projection matrices are a special class of matrices, closely related to the identity 
matrix, which may be singular. These matrices are useful when dealing with nonorthog- 
onal sets of vectors. 

Consider three linearly independent, nonorthogonal, unit column vectors ^[r/(i)^, 
i = 1,2, 3^ in a physical three space. Consider also three other column vectors ^ [?(!)) 
where i = 1,2,3^ which are related to the first three vectors by the relationship 


<?(D] 0?<j)> - Cy 


(Cl) 


where ( J denotes transpose of [[ / and 6y is the Kronecker delta defined by 


v° 

(i j) 

V 1 

(i = j) 


Next consider an arbitrary column vector [V) such that 

[v) = a 1 |r?(l)) + a 2 [t?( 2)) + a 3 |r?(3)) 


(C2) 


If equation (C2) is multiplied on the left by ^?(ij] , there follows from equation (Cl) 

H = (?«] [V> (C3) 

Equation (C3) gives a practical use for the jj(i)^ values; also the [?7(i))> and [j(i)^ 
terms are said to be mutually reciprocal. Next when [V) is written out in its expanded 
form, it is seen that 

[v> = ( fr(l)> <?(!)] + &(: 2)) (?(2)] + [r/(3)) (?(3)] ) [v> (C4) 

which shows that the identity matrix is defined by 

[i] = o?(i)) (?(i)] + &< 2 >> + &< 3 )> (?( 3 0 ( c5 ) 
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Let the term [rj( i)^ ^(ijj be dropped from [Y] and let Qv) be premultiplied 
by the result. It is seen from equations (Cl) and (C2) that [V^ remains unchanged 
except that the corresponding term a^ jrj(i)^ in equation (C2) is dropped. Thus when 
terms are dropped from [i], the resulting matrix is said to be the projection matrix for 
the directions corresponding to the terms retained in [Yj. 


The vectors [Y(i)^ may be found from the vectors jY(i)^> by forming a 

matrix [x] ^where the ith column of [x] corresponds to the components of |jj(i)^ 

and then forming the inverse of [x]. Then the components of the jth row of [x] 

are the components of |Y(j)^. Since j^?(i)^ terms are linearly independent, [X] 
exists. 


-1 

-1 
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